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We investigate the disordering of an initially phase-segregated binary alloy, due to a highly mobile 
defect which couples to an electric or gravitational field. Using both mean-field and Monte Carlo 
methods, we show that the late stages of this process exhibit dynamic scaling, characterized by a set 
of exponents and scaling functions. A new scaling variable emerges, associated with the field. While 
1-^ ' the scaling functions carry information about the field and the boundary conditions, the exponents 

^r^ , are universal. They can be computed analytically, in excellent agreement with simulation results. 






I. INTRODUCTION 



i-rt ', The kinetics of phase ordering in alloys, following a rapid temperature quench below the coexistence curve, has 

P3 ■ traditionally attracted much interest, both in the physics and material science communities 111. Remarkably, the 

O [ temporal growth of ordered domains, as measured, e.g., via a time-dependent structure factor, exhibits universal 

features, such as characteristic growth exponents and dynamic scaling, in analogy to critical phenomena. The "inverse" 

problem, i.e. the bulk disordering of a system after a rapid increase in temperature or in response to highly mobile 

defects, is also an interesting problem, occurring in the context of erosion or corrosion phenomena |0]. In addition, it 

is a testing ground for a fundamental question of statistical physics, namely, how a system approaches its final steady 

state, starting from an initial non- stationary condition. 

rvq ' In most real solids, microscopic atom-atom exchanges are mediated by defects, such as vacancies or interstitial sites 

(^ [p|. Thus, it is natural to describe such processes with the help of a three-state model, allowing for two species of atoms 

^~~^ and a certain number of empty (defect) sites [^ . Vacancy concentrations are typically extremely small (of the order 

j^ ' of 10~^). As a first step towards a deeper understanding of disordering dynamics, we will neglect any asymmetries 

— .^ ^ in the microscopies of the two atom species, leaving them for consideration at a later stage g. Thus, our study 

jrt ■ will be based on an Ising model with a very small admixture of empty sites, namely a single vacancy, undergoing 

Cj [ an ''^ upquencW from zero to infinite temperature. Since we wish to describe particles and holes, the dynamics will 

I ' conserve the number of each species separately. Only particle-hole exchanges will be permitted, in order to model the 

' O ' vacancy mechanism. We should note for completeness that a number of investigations have focused on downquenches 

^ •■ in similar models, i.e., vacancy- mediated domain growth [|6|. 

An alternate view of our study addresses the effect of a random walker on its background medium. Each step of 
the walker displaces part of the background, leaving tracks like a beach comber in the sand. These tracks can be 
monitored and display their own dynamics. In the simplest case M, the walker is a Brownian vacancy, exploring a 
lattice filled with particles of two species, labelled as "black" and "white", or "up" and "down" spins. Starting from 



5h ' a perfectly phase-segregated state, three distinct temporal regimes are observed, separated by two crossover times. 

. . 1 During the late stages of this process, the number of broken bonds exhibits dynamic scaling, characterized by a set of 
exponents and a scaling function. A mean-field theory allows for the analytic calculation of these features, in excellent 
agreement with simulations [Q. 

In order to test the range of universality of this model, we extend these studies to the case of a biased random 
walker: The vacancy or defect hops preferentially along one of the lattice directions, as if subject to gravity. In the 
language of electrostatics, the defect is charged, disordering an initial configuration of neutral atoms in the presence 
of an electric field. We emphasize that this is completely equivalent to having a neutral (or less massive) vacancy in 
a background of particles, all of which carry the same charge (or larger mass). As in the unbiased case, there is no 



feedback from the background to the defect, i.e., the motion of the vacancy is independent of the local background 
configuration. However, the choice of boundary conditions along the direction of the field now becomes important: 
we will consider two cases, one in which spatial inhomogeneities persist in the long-time limit, and another which 
approaches a homogeneous steady state. Intriguingly, in both cases the three temporal regimes noted before are 
still clearly observable. Our key goal is to explore to what extent scaling functions and exponents are universal, 
i.e., independent of boundary conditions and bias. Employing exact results, a mean- field theory and Monte Carlo 
simulations, we find that the scaling exponents are completely universal, but that the scaling functions can depend 
sensitively on the boundary conditions, through a new scaling variable which involves the bias. 

The paper is organized as follows. In Section II, we define our model, discuss the associated boundary conditions 
and introduce control parameters and observables of interest. We then turn to our findings. In Section III, we set 
the scene by characterizing the final steady states which are exactly known. On the basis of Monte Carlo data, we 
then demonstrate the emergence of three temporal regimes and the separation of time scales between them. This 
observation forms the basis of a mean-field theory, to be introduced in Section V. Its scaling predictions are compared 
to detailed Monte Carlo simulations in the following section. We conclude with some comments and open questions. 

II. THE MODEL. 

We begin by summarizing the discrete model underlying the Monte Carlo simulations. It is defined on a two- 
dimensional (2D) square lattice of dimensions L x L. Each lattice site is denoted by a pair of integers, r ={x, y), and 
can be occupied by either a black particle, a white particle or a vacancy. Multiple occupancy is forbidden. Thus the 
configurations of the model can be described by a set of spin variables {cr} which can take three values: a^ = +1 
(—1) for a black (white) particle, and a^ — for a vacancy. To model the minute concentrations of vacancies in real 
systems, we focus on the case of a single hole. The number of black {N^) and white {N^) particles is conserved and 
equal: N~^ = N~ = y K- ^^ ^^^ conclusions, we return briefly to the question of higher vacancy concentrations. 

The initial configuration is perfectly phase-segregated: Black and white particles each fill one half of the system, 
with a sharp flat interface between them, chosen to lie horizontally along the x-axis. The vacancy is located at the 
interface. Even though this structure is that of a ferromagnetically ordered state, the details of the particle-particle 
interactions are not important, since we will be interested in upquenches to infinite temperature, T = oo. Thus, the 
vacancy moves independently of its local environment. However, it is subject to a bias E, aligned transverse to the 
initial interface, i.e., upwards along the (positive) y-axis. This is of course equivalent to the particles experiencing 
an external gravitational or electric field pointing along —y. We stress that a factor 1/T has been absorbed into the 
parameter E; i.e., we are considering the high-temperature, high-field limit of a more complex interacting system. 

At each Monte Carlo time step (MCS), one of the four nearest neighbors of the vacancy is chosen at random. 
The exchange is performed using Metropolis [g| rates: min { 1 , e^"''^ } , modelling a /ocaZ electrostatic or gravitational 
potential. Here, Sy = 0, ±1 is the change in the y-coordinate of the vacancy, in units of the lattice constant a. Thus, 
moves against the field are exponentially suppressed while all others are accepted. No particle-particle exchanges are 
allowed. 

Next, we specify the boundary conditions. In combination with the bias, they determine whether any spatial 
inhomogeneities survive in the long-time limit, with potential consequences for scaling exponents or functions. In the 
unbiased case, the stationary state is completely uniform M , and the choice of boundary conditions affects at most 
nonuniversal amplitudes. To test the effect of spatial inhomogeneities, we consider two scenarios, namely, reflecting 
(also referred to as "brick wall", BWBC) or periodic (PBC) boundary conditions for the top and bottom edges. 
The right and left boundaries, being aligned with the drive, are not expected to play a significant role; we choose 
them to be periodic in both cases. These two scenarios differ in two important respects. First, reflecting boundary 
conditions allow some spatial inhomogeneities to persist, while periodic boundary conditions lead to homogeneous 
distributions. Second, we will see in the next section that, under BWBC, the system approaches an eguilihrium state 
in which all transport currents have subsided. In contrast, periodic boundary conditions are incompatible with a 
global Hamiltonian, so that the steady state is a non- equilibrium one, carrying a net mass current. 

The two control parameters for the Monte Carlo simulations will be the field strength, E, and the system size, L. 
To monitor the evolution of the system, we measure a "disorder parameter" , namely, the number of broken bonds 
(i.e., black- white nearest-neighbor pairs), A{L, E; t) m, as a function of time t. More detailed information is carried 
by the local hole and magnetization densities, defined respectively as 

(/.(r,i) = (5<,,,o) 

^(r,i) = (ar) (1) 



Here, (•) denotes the configurational average. The Kronecker-i5 ensures that lattice site r is the one occupied by 
the vacancy. Non-zero values of -0(1", t) indicate an excess of white or black particles at lattice site r, which is also 
a sensitive measure of the disordering process. The full time dependence of these densities can in general only be 
computed within a mean-field approach. However, their stationary forms are easily found from the known steady-state 
distributions, as we will presently see. 

We finally note that, for the analytic part of our work, it is straightforward to generalize our model to d dimensions: 
Denoting a lattice site by r ={xi, ...,Xd-i,Xd = y), the field selects a one- dimensional subspace, with reflecting or 
periodic boundary conditions, along the y-direction. The (d — l)-dimensional transverse subspace retains periodic 
boundaries. 

III. EXACT RESULTS: THE STEADY STATES. 

An exact solution of our model involves knowledge of the full time-dependent distributions, P{{ar}-,t), for the 
probability of finding configuration {iTr} at time t. This requires finding all eigenvalues and eigenvectors of the 
underlying master equation: 

a,P({ar},t) = ^ {W[{a'M'^r}]Pi{a'A,t) W[{ar}\{a',}]P{{a^},t)} . (2) 

IK} 

Here, VF[{iTr}|{o'r}] denotes the transition rate, per unit time, from configuration {cr^} into a new configuration {cr'^} 
which may differ only by one vacancy hop. In discrete time, it is just the Metropolis rate specified above. Even 
though Eqn (g) is only linear in the probabilities, a complete solution is usually feasible only for systems which are 
either very small {L < 3) or restricted to one dimension. Both cases are only of limited interest here. However, it is 
often possible to determine a particular eigenvector, namely the one associated with the (non-degenerate) eigenvalue 
zero: this provides us with the stationary limits, Podcr}) = linit^oo P({CTr}, i), of the full distributions. 

This procedure is particularly simple for the case of brick wall (reflecting) boundary conditions. Here, we can 
immediately write down the internal energy (the "Hamiltonian" ) of the system, being the electrostatic energy of 
a single charge in a uniform electric field. Thus, -Podcr}) (x exp ["^^ Eayd^^^o] is just the associated equilibrium 
Boltzmann factor. Since interparticle interactions are restricted to the excluded volume constraint, Pq is "color- 
blind" , i.e., it gives equal weight to all configurations of particles at fixed vacancy position. Of course, any additional 
interactions could easily be incorporated into the Hamiltonian. 

In contrast, the toroidal geometry of periodic boundary conditions, in combination with a uniform drive, prevents the 
existence of a global, time-independent potential. Therefore, the steady state is far from equilibrium. Yet fortunately, 
the stationary distribution is still exactly known pO| and even simpler: Pq is completely uniform, giving equal weight 
to all configurations of vacancy and particles. In this case, however, it is entirely unknown how to generalize this 
solution to more complex interparticle interactions. 

Given the stationary distributions, the steady-state densities, defined via 0o(r) = limi_,oo "/"(r, t) and 0o(r) = 
linif^oo 0^(1", i), are easily computed. For reflecting boundary conditions, (j)o is spatially inhomogeneous, following the 
usual exponential profile. Since the steady state is color-blind, the magnetization density is uniform. For simplicity, 
we place the coordinate origin into the center of the system. Thus, 

0o(r) = Cexp{Eay) and V'o(r) = , (3) 

where the normalization C = sinh(_Ea/2)/[L''~^ sinh(_Ea(L-|-l)/2)] ensures that the system contains only one vacancy. 
The second equation just reflects the fact that there are equal numbers of positive and negative particles |§]. It is 
noteworthy, and may not be entirely trivial at first sight, that the inhomogeneities in leave no trace at all in the 
magnetization distribution. 

The corresponding results for periodic boundary conditions are even simpler: 

0o(r) = j^ and tl^oir) = . (4) 

Here, (/)o(r) is also uniform, reflecting a single vacancy in a system of L'' sites. 

Finally, the saturation value for the number of broken bonds is readily found, since the particle configurations in 
the steady state are completely random. Since each of the four types of particle-particle bonds is equally likely, to 
leading order in 1/L, we obtain 



Asat{L,E)= \imA{L,E;t) = -L'' [1 + 0(1-^)] . 

The corrections include surface terms which distinguish the boundary conditions: for BWBC, we find 0{L^^) — 
—2'^^^/{dL) while the leading correction for PBC is only 1/i''. The vacancy plays an even smaller role since it affects 
only 2d bonds: this contribution is neglected here. 

Having established the key features of the steady state, we now turn to simulation results. 

IV. SEPARATION OF TIME SCALES. 

Our aim in this section is to motivate the key ingredient which will be required for most of our mean-field results, 
namely, that the vacancy equilibrates on a much faster time scale than the particles. Since this observation will be 
based on Monte Carlo data, we first summarize the parameters of our numerical work. 

We have performed detailed simulations of our model, using two-dimensional square lattices with L ranging from 
20 to 60. The lattice constant a is set to 1. Both types of boundary conditions, reflecting and periodic, have been 
implemented. For reflecting boundary conditions, the external field varies between and 1.0. Larger values of E or 
L equilibrate too slowly for our computational resources. For periodic boundary conditions, equilibration takes place 
much faster, and therefore, E varies between and 10.0. Our data are averaged over 10^ — 10^ realizations, depending 
on the desired accuracy. The length of the runs varies with system size, up to a maximum of 10® MCS. 






OMCS 



10^ MCS 



10^ MCS 






10^ MCS 



10*^ MCS 



10^ MCS 



FIG. 1. Sequence of snapshots showing the disordering process of 40a;40 system with E = 0.05 and BWBC. The black and 
gray squares represent the two types of particIes((T — ±1) and the white square denotes the vacancy(a — 0). The configurations 
were recorded after 0, 10^ 10"*, 10\ lO" and lO'^ MCS. 



To obtain a visual impression of the disordering process. Fig. 1 shows the evolution of a typical configuration, in a 
40 X 40 lattice with E = 0.05 and brickwall boundary conditions. The initial interface, completely smooth at t = 0, 
begins to break up as increasing numbers of particles are transported into the oppositely colored half of the system. 
Eventually, the system becomes completely disordered. Clearly, the number of broken bonds, A(L, E; t), is a suitable 
quantitative measure for the growing disorder, shown in Fig. 2 for the same set of system parameters. 
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FIG. 2. Plot of A{L,E;t) vs i, for L = 40, -E = 0.05 and BWBC. The inset shows the emergence of an early regime (I), an 
intermediate, or scaling, regime (II), and a late or saturation regime (III). See text for additional details. 

Similar to the unbiased case W, it clearly exhibits three regimes, drawn schematically in the inset: an early regime 
(I), an intermediate regime (II), and a late or saturation regime (III) in which the system has reached steady state. 
As the system size increases, the intermediate regime spans a widening time range, suggesting that the three regimes 
are temporally well separated. In practice, this is already the case for L > 20. This separation of time scales is 
also observed for larger values of the bias. Physically, the three regimes are characterized as follows. For early times 
(regime I), the vacancy is still localized in the vicinity of its starting point. Given the properties of random walks, 
this regime is expected to depend strongly on the dimensionality. After a time of 0{L^), however, the vacancy has 
explored the whole system and is effectively equilibrated. This marks the onset of the intermediate regime. In contrast, 
the particle distribution is still strongly inhomogeneous, equilibrating only upon entering the saturation regime. A 
completely analogous picture emerges for periodic boundary conditions. 

Our goal in the following is to test for dynamic scaling. Since our data involve the three variables L, E and t, 
as well as two types of boundary conditions, the appropriate scaling forms may not be a priori obvious. We shall 
therefore first consider an analytic approach, before turning to a detailed comparison with the data. 



V. MEAN-FIELD THEORY 



A. The equations of motion. 



Since an exact solution of the full master equation is not available, we seek a simpler formulation of the dynamics. 
For our purposes, a set of equations of motion for the two conserved densities, (j) and ip, provides a suitable starting 
point. We immediately focus on the d-dimensional case. A spatially discrete version of the desired equations is easily 
obtained from the master equation, via the definition (m). Not surprisingly, the excluded volume interaction mixes 
correlation functions of different orders, generating the usual BBGKY hierarchy. To obtain a closed set of equations, 
a mean-field approximation is required: higher correlations are approximated by products of appropriate single-point 
averages. Finally, we take the continuum limit, by letting the lattice constant a vanish at fixed physical system size 
L and field E. The microscopic time scale is identified with a^/2d. Thus, the resulting densities are functions of a 
continuous d-dimensional coordinate r = {xi, . . . , x^-i, Xd = y) and time t. The spatial origin is chosen at the center 
of the system so that — L/2 < a;^ < L/2 for i = 1,2, ...d. Letting y denote the unit vector along the direction of E, 
we finally obtain the basic equations of motion: 



dt(l){v,t) = V[V^{v,t) - Ey(l){v,t)] 

dt^ir, t) = V[0(r, t)V7/-(r, t) - ^(r, t)V0(r, t) + Ey^{v, t)^{y, t)] 



(5) 
(6) 



Here, V is the d-dimensional gradient. Clearly, both equations take the form of continuity equations, reflecting 
the conservation laws. The right-hand sides are the (negative) gradients of the associated hole and magnetization 
currents. In the field- free limit E = 0, these are just the results of Ref. ITJ. Of course, the bias induces additional 



terms which reflect systematic Ohmic currents: Since particles can only move when the vacancy is present, the extra 
magnetization current, —Ey(j)ip, is proportional to (j). In contrast, the hole current Ey(f> is independent of the local 
particle background, as one should have expected. We also note here that E has dimensions of inverse length. This 
suggests that the diniensionless combination LE will emerge as a convenient scaling variable. 

Eqns (ph and (o) have to be supplemented with the constraints on the total hole and particle ||] numbers: 



(r, t) = 1 and 



V'(r,f) =0. 



(7) 



Here, V — L'^ \s the volume of the system. Next, we consider the boundary conditions. In the PBC case, the solutions 
of Eqns (g) and (||) must be fully periodic functions of r, with period L. For reflecting boundary conditions, the 
periodicity of the solutions is restricted to the transverse subspace. Along E, we demand instead that the hole and 
magnetization currents vanish on the boundary y = ±L/2. 

Finally, we specify the initial conditions. The vacancy starts at the origin 

0(r,O)=,5(r) (8) 

and the particles are perfectly phase-separated: 

V^(r,0)==2%)-1 (9) 

This completes the discussion of the equations of motion and their constraints. We now focus on their solution. 

B. Solutions. 

As a starting point, we first establish the steady-state solutions of Eqns (|^) and (H). When turning to the full time 
dependence, progress is only possible if we simplify the equations of motion by invoking the separation of time scales. 
We shall mostly focus on brick wall boundary conditions, since this situation is physically more complex, as we shall 
shortly see. At the end, we briefly summarize our results for the fully periodic case. 




FIG. 3. Simulation data for the steady state hole profile with L = 20 and £ = 0.01(+), 0.05(x), 0.10(*) and 0.15(D), BWBC. 
The dashed lines denote the mean-field profiles, EqndlO) and (pd[), at the same parameter values. 

Anticipating inhomogeneities along the w-direction only, due to the bias, we seek steady-state solutions of the form 
(/)(r, t) = 4>o{y) and ip{r, t) — ipo{y). Eqn (|5|) can now be integrated once, and the reflecting boundary conditions force 
the integration constant, i.e., the hole current, to zero. We can immediately integrate again, whence 



0o(y) =Cexp(^2/) 
The constant C follows from the normalization condition, ([7|): 

EL 



C 



2L''smh{EL/2) 



(10) 



(11) 



Not surprisingly, this expression is just the continuum limit of the exact result, Eqn (^. One should note the emergence 
of the conjectured scaling variable LE. In Fig. 3, we compare the mean-field profiles for L = 20 and a range of E's 
with the corresponding Monte Carlo data. The agreement is clearly excellent. Finally, a similar integration of Eqn 
(^, using (|lO|), yields ipoiv) = 0, also consistent with the exact form. 

Unfortunately, the solution of the time-dependent equations is not so simple. The full solution 0(r, t) of Eqn (H) 
must be found and inserted into Eqn (ph. The latter poses a formidable problem: even though it is a linear, second 
order partial differential equation of parabolic type, its coefficients are rather complicated functions of space and time. 
Fortunately, a complete solution is possible if we restrict our attention to the intermediate and late regimes. Here, 
the vacancy density has already reached steady state, so that only the purely exponential (poiv) enters into Eqn m). 
Of course, the reduced equation should be supplemented with the initial condition (O). Symmetry considerations as 
well as the data suggest that the magnetization density is uniform transverse to E^ so that we seek a solution of the 
form ip{r,t) — ip{y,t). Letting d denote the partial derivative with respect to y, we arrive at a much simpler version 
of Eqn dg) , describing the magnetization density at intermediate and late times: 



It is subject to the initial condition 



and the brick wall boundary condition 



d^i±^,t) 







(12) 



(13) 



which ensures that the magnetization current through the ends of the system vanishes. 

This differential equation can be reduced to an eigenvalue problem with a complete and orthonormal set of eigen- 
functions {Un{y)} which are linear combinations of Bessel functions. The eigenvalues k„ are real, discrete and strictly 
positive. Deferring all mathematical details to the Appendix, we only quote the form of the solution, expressed as an 
eigenfunction expansion: 



Hy,t) 



^ AiUniv) exp(-K„i) 

n=l 



(14) 



The expansion coefficients An are chosen so as to match the initial condition (S). 




FIG. 4. Simulation data for charge density profile with L = 20, E = 0.1 and BWBC after 10'*(+), 5a;10''(x), 10^ {*\3xl0^ (a) , 
5a;10''(H), 7xl0^(o), and 10^(«) MCS. The dashed lines show the corresponding mean-field profiles based on Eqn(|l4). See text 
for details. 



Combined with the explicit expressions for its constituents (cf. the Appendix), Eqn (14) is exact within our mean- field 
approach. It agrees remarkably well with Monte Carlo data, as demonstrated by Fig. 4. There, we compare the sum 
of the first 1000 terms in Eqn (14), for several values of t, with simulation results. A single fit parameter is required, 
which converts t into the number of MCS. One sees clearly how the initially phase-segregated system disorders with 
time. 

Having established the structure of the densities, we turn to the disorder parameter. Being the number of broken 
bonds, it is directly related to the (nearest-neighbor) Ising "energy" of our model, (~J2<rr'>'^r(^r')- Since the 
continuum limit of the latter is proportional to Jy ip'^ ||ll[ , we obtain Q : 



7 



A{L,E-t)^-L'' 



1 



-L/2 



L/2 



dyip^iy,t) 



+ OiL''-') . 



(15) 



We emphasize that the boundary conditions only affect surface terms. These will be neglected in the following, leading 
only to small errors since L > 20 in the data. To the same accuracy, the initial condition for A is just A{L,E;0) 
= 0. Using Eqn (Um and carrying out the integral, we find the time evolution of the disorder parameter, for brick 
wall boundary conditions: 



A„{L,E-t) = -L'' 



-. oo 

1 - -— ^ ^2 exp(-2K„i) 



n=l 



(16) 



Here, we have introduced a subscript, R ("reflecting") or P ("periodic"), to distinguish the two types of boundary 
conditions. 

We conclude by focusing briefly on the case of fully ■periodic boundary conditions. The steady-state solutions of 
Eqns (||) and (|6|), consistent with the constraints (|^), are </'o(r) = 1/L'^ and V'o(r) = 0. To obtain the time-dependent 
excess density 'ijj{r,t), we invoke the separation of time scales again. Noting that the vacancy density has equilibrated 
at the beginning of the intermediate regime, we insert its steady-state form, (j)o, into Eqn (0), whence 

5tV'(r,t)=V[0o(V+£;y)^(r,i)] 

This simplified form holds in the intermediate and late regimes. A Galilei-transformation r = r' + E(f)oyt recasts it 
as an ordinary diffusion equation, with diffusion coefficient ^c,: 



at^(r',t) = 0„V2V'(r',i) 



(17) 



Quite remarkably, we observe that all effects of the field disappear in a suitably chosen co-moving frame. Thus, the 
solution of Eqn (O) and the accompanying disorder parameter can be read off immediately u^ from Ref. W. In 
the original frame, the whole profile ip{r,t) drifts with velocity —EL~'^y. Of course, no such drift is observed in 
Ap{L, E; t), by virtue of the spatial integration in Eqn dlq): 



ApmE;t) = -L'' 



- y 

n=l,3,. 



exp(— 2K„i) 



(18) 



where k„ = {2TTnf / L''+^ . 

This completes the discussion of our mean-field equations and their solutions. We now turn to the underlying 
scaling behavior and compare it to Monte Carlo data. 

VI. DYNAMIC SCALING: ANALYTIC PREDICTIONS AND NUMERICAL TESTS. 



In this section, we extract the scaling forms for the number of broken bonds, A{L, E; t). We begin with Eqn (Hq) 
for brick wall boundary conditions. In the Appendix, it is shown that both the coefficients An and the eigenvalues k„ 
exhibit characteristic scaling forms: A^ is a function of the new scaling variable LE alone, while the eigenvalues obey 

«;„ = L-('^+2)g„(Li?) . 

Therefore, the scaling form of Ar{L, E; t) is apparent, namely, 

AR{L,E-t) = L''TR{LE,t/U). (19) 

where J-r is an appropriate scaling function. The temporal scale factor t^ is itself still a function of LE and can be 
chosen in different ways. Here, we focus on the crossover from the intermediate to the saturation regime and define 
tn as a measure for the late crossover time. 



t,{L,E) = K^'^L''+'TciLE) 



(20) 



The limits of the scaling function Tc = gi are discussed in the Appendix: It approaches its field- free limit for vanishing 
LE and increases exponentially with LE in the opposite limit, reflecting the increasing confinement of the vacancy. 



We already note one of our key results here, namely the emergence of a new scaling variable, LE. Its physical origin 
is clear: it determines how easily the vacancy can escape from the top {y = L/2) edge of the system, where the field 
tends to localize it, in order to travel to the center of the system where most of the disordering is taking place. 

It is interesting to contrast the behavior of the disorder parameter, and hence the scaling function Tji for "late" 
and "early" times, t/tc ^ 1 and t/tc <C 1 (but already within the intermediate regime) respectively. For late times, 
the disorder parameter has saturated so that liuiy^oo ^R{x,y) — d/2 independent of a; = LE. The short-time limit 
requires some care since Eqn dlq) does not converge well there. However, using a Poisson resummation (see the 
Appendix for details), we can recast it in an alternate form that converges rapidly for small C oc t/tc'- 



^H(L,i?;i) ~ ^L'^^V^Il + O [cxp(-7rVC')] } 



(21) 



Thus, we conclude that the disorder parameter increases as a power law, t^ , for short times, with an exponent /3 = 1/2. 
The short-time scale factor exhibits a similar scaling form as the late crossover time. 



t,{L,E) = L^+'T,{LE) 



(22) 



but we emphasize that its scaling function Ts differs from Tc in its dependence on LE. This feature finds its origin 
in the short-time limit of the scaling function, \m\y^Q!FR{x,y) = GR{x)y^''^, expressed by Eqn (pi|): The presence 
of the second argument, x = LE, leads to the non-trivial prefactor Gr{x) which translates between the two scaling 
functions t^ and Tc. 

Before turning to a comparison with the data, we first contrast Eqns (|l9|-P2|) with the corresponding results for 
periodic boundary conditions: 



ApiL,E;t) ^L'^TpitL 



-(d+2) 



L-'Tpit/t,). 



(23) 



All dependence on E disappears here: The scaling function Tp can be found in Ref. |7|, and the late crossover time 
tc is just L'^+^. A Poisson resummation yields the short-time behavior Ap oc ^L'^^Jt|tc{l + O [exp(— tt^/C^)] }, for 
C, oc t/tc "^ 1- Thus, we also find a power law increase here, Ap oc t^ , with the same exponent /3 — 1/2. 
These predictions are tested in Figs. 5 and 6, for periodic and reflecting boundary conditions, respectively. 
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FIG. 5. Plot of the total number of broken bonds with PBC. (a) Plot of A{L,E\t) vs t of 40x40 with 
E = 0.00(+),0.10(x), 1.00(*) and lO.O(n). (b) Scaling plot of A{L,E-t)/L'^ vs t/L'^ for L = 30 - 60 with E = 0.00, 1.00 and 
10.0. 



In Fig. 5a, we show the disorder parameter Ap, for one system size but several values of E. Some slight deviations are 
observable at short times, before the vacancy density equilibrates. In the intermediate and late regimes, however, it 
is quite striking, but entirely consistent with our expectations, that all data points fall onto the same curve, without 
any scaling being required. In Fig. 5b, we show a scaling plot for several L and E. Excellent data collapse is observed 
in the intermediate and late regimes, where our scaling predictions are expected to hold. Again, we emphasize that 
the scaled axes depend only on L, but not on E. 
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FIG. 6. Plot of the total number of broken bonds with BWBC. (a) Plot of A{L,E\t) vs 
E = 0.00(+),0.10(x),0.20(*), 0.30(D) and 0.50(B). (b) The scaling plot of A{L,E-t)/L'^ vs t/L" for L= 
and 60(a) with LE = 2.0. 
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Fig. 6a shows the disorder parameter Ar for reflecting boundary conditions. Several values of E are plotted, for 
L — 40. As anticipated, the mixing process proceeds more slowly for larger values of E, since these tend to confine 
the vacancy more strongly to the upper edge (y = +L/2) of the system. The corresponding scaling plot is shown in 
Fig. 6b: Here, several values of L and E are chosen such that LE = 2.0 remains constant. The data collapse in the 
intermediate and late regimes lends full support to our prediction, Eqn (^9). 

Focusing on the intermediate regime {0{L^) < t <C tc), Figs. 5b and 6b show very clearly that the disorder 
parameter increases as a power law there, A{L,E;t) oc t^ . Since this behavior persists over at least three decades, 
we can extract a reliable numerical estimate for the exponent, resulting in /? = 0.5 ± 0.01 for both reflecting and for 
periodic boundary conditions, in agreement with our prediction. 

Our results so far can be summarized more succinctly. Defining a set of characteristic exponents and scaling 
functions, the scaling of the disorder parameter, for both types of boundary conditions, can by written as 



AiL,E;t)r^ L°'T{LE,t/t^) with t^{L,E) ^ L^t{LE) 



(24) 



Remarkably, we find that the exponents a and z are universal^ i.e., independent of boundary conditions and drive. 
Our results, a — d which is exact and z — rf + 2 which follows from our mean-field theory, are completely consistent 
with the Monte Carlo data and agree with the corresponding exponents for the unbiased case 0. The short-time 
behavior can be written as 



AiL,E;t)^L''{t/T,) 



P 



(25) 



with a "growth" exponent j3 = 1/2 which is also manifestly universal. The additional exponent a is related to the 
others via the consistency condition A{L, E; tc) ~ A{L, E;t ^ oo), whence a = d— zf] = [d— 2)/2. This scaling law 
is obviously satisfied by our results. 

In contrast to the exponents, the scaling Junctions T and t are profoundly affected by the bias and the boundary 
conditions. For PBC, we just recover the results of the unbiased case, whereas a nontrivial dependence on a new 
scaling variable, L_E, emerges in the brick wall case. 

VII. CONCLUSIONS 

In this paper we have analyzed the vacancy-mediated disordering of a binary alloy, in response to an upqucnch 
from zero to infinite temperature. The system is placed in a gravitational or electric field, and two types of boundary 
conditions, reflecting ("brick wall") and periodic, are studied. Starting from a perfectly phase-segregated initial 
configuration, the vacancy mediates atom exchanges, leading to fully disordered exactly known final states. For 
brickwall boundary conditions, the final state is an equilibrium one, characterized by a Boltzmann distribution and 
an exponential hole density profile. In contrast, the steady state established by periodic boundary conditions is a 
non- equilibrium one, with homogeneous configurational distribution and profiles. Using Monte Carlo simulations, we 
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monitor the time evolution of the system by measuring local profiles and a disorder parameter, i.e., the number of 
broken bonds. For both types of boundary conditions, three temporal regimes are identified: an early one, where 
the vacancy has not reached the boundaries of the system, an intermediate one in which the vacancy has already 
equilibrated, but the particle distributions are still inhomogeneous, and a late regime where the system has reached 
steady state. To predict scaling exponents and scaling functions, we develop a theoretical description in terms of a 
set of mean-field equations of motion for the local densities. Invoking the separation of time scales, the mean-field 
equations can be solved exactly in the intermediate and late time regimes, providing us with explicit expressions 
for the observables of interest. We find that the disorder parameter exhibits dynamic scaling and observe excellent 
agreement of Monte Carlo data and mean- field predictions. Our key result is a set of universal scaling exponents, 
independent of bias and boundary conditions. For example, in all cases the number of broken bonds increases as i^'^, 
before saturating at a value of 0{L'^). In contrast, the scaling functions for the disorder parameter exhibit a lesser 
degree of universality: While they remain independent of the bias in the case of periodic boundary conditions, an 
additional scaling variable, LE, must be taken into account for reflecting boundary conditions. This variable also 
controls the shape of the excess density (magnetization) profile. 

So far, we have restricted our attention to systems containing only a single vacancy. However, in real systems, one 
should expect that the number of vacancies, M, scales with system size [[7|. For generality, we introduce the vacancy 
number exponent 7, with < 7 < d, so that M ~ L'^ . Thus, 7 = d — 1 describes a situation where the vacancies 
("defects") prefer the interfacial region before the upqucnch occurs. Our previous results correspond to the case 7 = 
but are easily extended. We simply need to modify the normalization condition for the hole density profile, Eqn (|^), 



so that the constant C, given in Eqn (11), becomes C — EL/[2L'^ ' sm'h{EL/2)]. As a consequence, the late crossover 
time, Eqn (po|), now scales as tc{L,E) = L'^^^^'^rdLE), so that z = d + 2 — 7 and a = (d + 7 — 2)/2. Of course, these 
exponents are again universal. 

We conclude by noting the different symmetries which characterize the BWBC and the PBC case, in the presence 
of the bias. Periodic boundary conditions are compatible with translation invariance, but violate the detailed balance 
condition: the driving force is not compensated by a chemical potential gradient. In contrast, detailed balance 
holds for reflecting boundary conditions, but translation invariance is broken. It is quite remarkable that the scaling 
exponents remain unaffected by such profound differences in symmetry. 
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APPENDIX 

A. Full solution for the magnetization density. 

In this appendix, we present the mathematical details of solving Eqn (|l2| ) for brickwall boundary conditions, 
restricted to the intermediate and late stages of the disordering process, i.e., 0{L^) <C t. Thus, in the following, 
"early" refers to times at the onset of the intermediate regime while "late" times are deeply within the saturation 
region. To recall, we seek a solution, ^p{y,t), to the partial differential equation 

dMy,t) = d[My)d^iy,t)] (26) 

where 

, , EL . . 

'^°(^)=2L'^sinh(i?i/2)^"P(^^) 

subject to the initial and boundary condition V'(y)O) — 26{y) — 1 and dip{±^,t) — 0. First, we introduce the new 
variable 

y (Po[y') E(j)o{y) 



which is strictly positive. This reduces Eqn (|2g) to a diffusion equation with a spatially varying diffusion coefficient: 

dt^(x,t)^Exdli^{x,t). (28) 



subject to the boundary condition dx'4'(x±,t) = at a;± = [E(J3o{^L/2)] . We note that the definition (|2^) implies 
< x+ < X-. Next, we separate variables, according to the ansatz 

^:{x,t)^T{t)f{x). (29) 

whence we obtain two ordinary differential equations: 

7rTn 

-— + nT = Q (30) 

at 



dx'^ Ex 



2 ■ n.^/ = 0' (31) 



The constant k must be positive in order to be consistent with the steady state solution, lim^^oo i^iy, t) = i^o{y) = 0. 
The first of these equations describes a simple exponential decay. The second constitutes a well-defined Hermitcan 
eigenvalue problem, with eigenvalues k and eigenfunctions /. For convenience, we define k/ E = a^/4 and transform 
( pf| ) into the differential equation for the Bessel functions [^ , via u = a^/x. The solutions are 

T{t) = exp(-Ki) (32) 

f{x) — Aa^/xJi{a^) + Ba\/xNi{a\/x). (33) 

Here, A and B are integration constants. We have two boundary conditions, one at each end of the system. Using 
the recursion relations for the Bessel functions, we can eliminate one of the integration constants, e.g., 

B = -A 2\ ^ ', 34 

and specify the allowed eigenvalues a„ as the solutions of the implicit equation 

= Nii{a.a^fxZ)Ji^{an^fx^) - A^o(c^n^/i^+)-^o(an\/iE^) (35) 

= iVo(Az„) Jo(z„) - iVo(^n) Jo(Az„) (36) 

The second line, with Zn = otnJxj^ and A = 4/^^ = exp(Lii^/2) > 1, is a more standard form of the eigenvalue 
equation [O. Since both Jq and Nq oscillate, this equation has infinitely many solutions. The eigenvalues are real, 
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non-degenerate and discrete; they increase monotonically with n. The lowest ones are easily determined numerically 
for different A. For large mr/{X — 1), there is an asymptotic expansion, 



nvr I 1 /A-1 

A- 1 I ^ 8A V nn 



O 



A-1 



(37) 



To find the cigcnfunctions {[/„(a;)}, we need to normalize the /'s. For convenience, we introduce the auxiliary function 

F{s) = s[No{any/^)Ji{s) - Jo(a„V^)iVi(s)] (38) 

and define the normalization constants 

c = ^ ^ ^ (39) 

Then, the eigenfunctions take the form: 

Unix) = -^F{anVx) 
an 

= CnVx [No{an^/xZ)Ji{an^/x) ~ Jo{an^/xZ)Ni{an^/x)] . (40) 

They are real and orthonormal: 

Snm = / dxx^^Un(x)Um{x) ^ E I dy [/„(?/) C/„(y) 

Jx+ J-L/2 

where the second equality expresses the orthonormality condition in terms of the original variable y. They form a 
complete set so that the full solution for the magnetization density can be written as an expansion: 



V'(a;,t) = ^A„t/„(x)exp(-K„i) , 



(41) 



where k„ = i?a^/4, and x — [Ecjyoiy)] ■ The expansion coefficients A^ are chosen such that the initial condition is 
satisfied. With Xo = [Eipoi^)] , this yields: 



4c 
An = —^ [NQ{an^/x^)Jo{an^/x~) - N^){an^/x~)Jo{an^/x^)\ 



A^o(^/Az„) Jo(Az„) - iVo(Az„) Jo(\/Az„) 



[F2(Az„)-F2(z„)]i/' 



(42) 
(43) 



This completes the solution. Of course, it is given rather implicitly in terms of the eigenvalues. To make our 
expressions more transparent, we track the key system parameters, L and E^ through these manipulations, in order 
to exhibit the scaling properties of the theory. 

B. Scaling analysis. 

1. Eigenvalues and expansion coefficients. 

We first establish the scaling properties of the eigenvalues, k„. Beginning with Eqn (pq), we conclude that the z„'s 
are functions of A alone. Since A — exp{LE/2), each Zn depends only on the scaling parameter LE. To obtain a 
similar conclusion for k„ = Eaf^/A = £'z^/(4a;+), we recall that 



l/x+ = EM+L/2) = 



Thus, the desired scaling form for the eigenvalues is 



(EL)^ exp{EL/2) 
2Ld+i smh{EL/2) 
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«„ == L-'^'^+^^g^LE) (44) 

where the scalmg function gn{x), with x = LE, is given by 

g„(x)EEz^(.)^!^fi^ (45) 

8smh(a;/2) 



Its hniits for small and large argument are easily found from Eqn (36) and the asymptotic form ( 



lim (7„(x) = (mr) [I + 0{x)] and lim gn{x) ~ a;„ exp(— x/2) , (46) 

where x„ denotes the n-th zero of the Bessel function Jo- 
in particular, we are interested in the late crossover time, Eqn {Uw, defined as the inverse of the first eigenvalue, 

t,{L, E) = K^^ = L''+\,{LE) . (47) 

This characteristic time is a measure for when the crossover from the intermediate to the saturation regime occurs. The 
behavior of its scaling function, Tc{x) = g^ (x), follows immediately from Eqn ([46[): \imx-tQTc{x) = 7r~^ [1 + 0{x)] 
and \imx~^oo Tc{x) ~ (0.1729.../a;^) exp(a;). 

Finally, we will need the scaling behavior of the expansion coefficients An. From Eqn ([l3[), it is immediately 
apparent that these coefficients depend only on the combination LE so that An — An{LE). 

2. The Poisson resummation of the disorder parameter: short-time behavior. 

While Eqn (Hq) for A{L,E;t) is of course completely exact within our mean-field theory, it converges rapidly only 
for late times, K„i 3> 1. There, keeping only the n = 1 term in the sum already results in an excellent approximation. 
In contrast, Eqn dlq) is not very practical if we wish to extract the observed power law at early times. Fortunately, 



a Poisson resummation \ 15 of Eqn (|l6| ) allows us to recast the disorder parameter in a form that converges rapidly 



in the short-time limit. Some details of this procedure form the content of this section. 



We recall Eqn l\i§ 



AR{L,E-t) = '^L' 



^ OO 

1 - -— ^ Al exp(-2K„t) 



(48) 



where both An and k„ depend on the summation index via the eigenvalues z„. The key to the Poisson resummation 
resides in the following three statements: 

First, from the discussion below Eqn (p^, we recall the initial condition on the disorder parameter, namely 
Ar{L,E-{)) = 0. This implies 



^EA^^-1' (49) 



EL 

Tl=l 

so that there is no constant term in the short-time expansion of Ar. 

Second, considering any finite number of terms in Eqn (Es) can only generate a /mear time dependence, Ar{L, E; t) oc t. 

Therefore, the anticipated short-time behavior Ar, cx y/t must be controlled by the large n contributions to the sum. 

Hence, these are crucial for our purposes. 

Third, for sufficiently large n > rio, we can always approximate the eigenvalues Zn by their explicit asymptotic form, 

Eqn (|37|). Since the latter holds provided mr/[X — 1] '3> 1, the critical no ~ A = exp{LE/2) increases rapidly with 

LE. However, this presents no problem since only finite values of LE are of interest to us. 

In summary, to obtain the short-time behavior of Ar it is sufficient to replace the eigenvalues z„ by their asymptotic 
form everywhere. Any errors generated in this way are at most linear in t. In this manner, the dependence on the 
summation index n becomes explicit, and we can apply the Poisson resummation formula Pq |, which holds for any 
continuous, bounded function f{x), provided its Fourier transform F{uj) = 2 J dxf{x) cos{ojx) is well-defined: 



2/(0) + E /("O = - 2^(^) + E ^(^) 

n— 1 L m— 1 ) 



(50) 
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Here, <; is the parameter that controls the convergence of the sums. In our case, we identify 



'-Y(A-1)3(A + 1) 
The resummation is now straightforward and results in 



^ ' L-id+2)t (51) 



An{L, E- 1) = ^i-^-i^L^c {1 + O [eM--'ie)] } 



2 v^ 
where the scale factor for the short-time scaling is given by 



,<i+2 4sinh(Li?/2) 



tJL, E) = 1"+' \ — Vt- (52) 

^ ' ' LEcxp{LE/2) ^ ' 



While this characteristic scale obeys the same scaling form as the late crossover time, 

U{L,E) = U^-+\,{LE) with r..(x)^ ;^"Pf//) 

4smh(.T/2) 

we note that the scaling function Ts{x) is different from Tc{x). 



[1] J.D. Gunton, M. San Miguel and P.S. Sahni, in: Phase Transitions and Critical Phenomena, Vol. 8, eds. C. Domb and 

J.L. Lebowitz (Academic, New York, 1983); H. Furukawa, Adv. Phys. 34, 703 (1985); K. Binder, Rep. Prog. Phys. 50, 

783 (1987), and in: Phase Transition of Materials, eds. R.W. Cahn, P. Haasen and E.J. Kramer (VCH, Weinheim, 1991); 

J. S. Langer, in Solids far from Equilibrium, ed. C. Godreche (Cambridge University Press, 1992); A. J. Bray, Adv. Phys. 

43, 357 (1994). 
[2] H-F. Meng and E.D.G. Cohen, Phys. Rev. E51, 3417 (1995); J. Krug and P. Meakin, Phys. Rev. Lett. 66, 703 (1991); 

R. Cuerno, H.A. Makse, S. Tomassone, S. Harrington and H.E. Stanley, Phys. Rev. Lett. 75, 4464 (1995) and references 

therein; S. G. Corcoran, to appear in: Proceedings of the Electrochemical Society, Symposium on Critical Factors in 

Localized Corrosion III, 1998. 
[3] see, e.g., W.D. Callister Jr., Material Science and Engineering: An Introduction (Wiley, 1994). 
[4] M. Blume, V.J. Emery and R.B. Griffiths, Phys. Rev. A4, 1071 (1971); R.B. Potts, Proc. Camb. Phil. Soc. 48, 106 (1952); 

F.Y. Wu, Rev. Mod. Phys. 54, 235 (1982). 
[5] B. Schmittmann, W. Triampo and R.K.P. Zia, to be published in Brazilian Journal of Physics (2000). 
[6] O.G. Mouritsen and P.J. Shah, Phys. Rev. B40, 11445 (1989); P.J. Shah and O.G. Mouritsen, Phys. Rev. B41, 7003 

(1990); K. Yaldram and K. Binder, Z. Phys. B82, 405 (1991) and references therein; K. Yaldram, G.K. Khahl and A. 

Sadiq, Solid State Comm. 87, 1045 (1993); C. Frontera, E. Vives, T. Castan and A. Planes, Phys. Rev. B51, 11369 (1995) 

and references therein; S. Puri, Phys. Rev. E55, 1752 (1997); P. Fratzl and O. Penrose, Phys. Rev. B50, 3477 (1994) and 

Phys. Rev. B55, R6101 (1997); S. Puri and R. Sharma, Phys. Rev. E57, 1873 (1998). 
[7] Z. Toroczkai, G. Korniss, B. Schmittmann and R. K. P. Zia, Europhys Lett. 40, 281 (1997). 
[8] This equation is strictly correct only for systems with an odd number of sites. For a system with an even number of sites, 

the exact relation is A''+ = N~ + 1 = -y, i.e., iV+ — N^ + 0(L^^). However, this correction will be neglected in the 

following: even for our smallest system (20 x 20), it is below 1%. 
[9] N. Metropolis, A.W. Rosenbluth, M.M. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953). 
[10] F. Spitzer, Adv. Math. 5, 246 (1970). 

[11] D.J. Amit, Field Theory, the Renormalization Group, and Critical Phenomena (World Scientific, 1984). 
[12] Ref. [H uses a slightly different normalization. 
[13] see, e.g.. Handbook of Mathematical Functions, with Formulas, Craphs, and Mathematical Tables, eds. M. Abramowitz and 

I. A. Stegun, (Dover, New York 1974). 
[14] To the best of our knowledge, this is a new sum rule for Bessel functions. 
[15] D.S. Jones, The Theory of Generalised Functions (Cambridge University Press, 1982). 



15 



